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Abstract 

Methanol and its precursor formaldehyde are among the most studied organic molecules in the interstellar medium 
and are abundant in the gaseous and solid phases. We recently developed a model to simulate CO hydrogenation via 
H atoms on interstellar ice surfaces, the most important interstellar route to ILCO and CH3OH, under laboratory 
conditions. We extend this model to simulate the formation of both organic species under interstellar conditions, 
including freeze-out from the gas and hydrogenation on surfaces. Our aim is to compare calculated abundance ratios 
with observed values and with the results of prior models. Our model utilises the continuous-time, random-walk 
Monte Carlo method, which — unlike other approaches — is able to simulate microscopic grain-surface chemistry 
over the long timescales in interstellar space, including the layering of ices during freeze-out. Simulations under 
different conditions, including density and temperature, have been performed. We find that H2CO and CH3OH form 
efficiently in cold dense cores or the cold outer envelopes of young stellar objects. The grain mantle is found to have a 
layered structure with CH3OH on top. The species CO and H2CO are found to exist predominantly in the lower layers 
of ice mantles where they are not available for hydrogenation at late times. This finding is in contrast with previous 
gas-grain models, which do not take into account the layering of the ice. Some of our results can be reproduced by 
a simple quasi-steady-state analytical model that focuses on the outer layer. Observational solid H2CO/CH3OH and 
CO/CH3OH abundance ratios in the outer envelopes of an assortment of young stellar objects agree reasonably well 
with our model results, which also suggest that the large range in CH 3 OH/H 2 observed abundance ratios is due to 
variations in the evolutionary stages. Finally, we conclude that the limited chemical network used here for surface 
reactions apparently does not alter the overall conclusions. 



1 Introduction 



Methanol and its precursor formaldehyde are among the most studied organic molecules in the interstellar medium 
(ISM). Both species have been observed in the gas phase and as constituents of ice mantles. Whereas gas-phase detec- 

lj0rgensenetaill2OO5t iMaret et alJl2004l2005tlSchoier et al.ll2 004: 



tions of both molecules have been numerous, (e.g., 

van der Tak et al. 2000) and methanol is a well-known constituent of the ice lAH amandola et al 



1992t IPartois et al 



1999; Her bst & van Disho eck 2009), there are only a few secure solid-state detections of H2CO. These detections 
occur mos tly in high-mass young stellar source s such as W33 A, GL2136, and GL 989 , and often only upper limits 
are given (Boogert et al. 2008; Gibb et al. 2004 ). For example. IPontoppi dan et al. (2004) found an upper limit for the 
H2CO ice abundance in the outer envelope of the low-mass young stellar source Serpens SMM 4 of 5% with respect 
to water ice. CH3OH ice has a much higher abundance of 28% with respect to water in this source, resulting in an 
H2CO/CH3OH ice ratio of less than 0.18. Typically, in sources where both species have been detected, the observed 
H2CO/CH3OH ratio is still smaller than unity, with values ranging from 0.09 to 0.5. Sub-millimetre observations of 
the gas in seven massive hot-cores show a similar trend: iBisschop et al] d2007l) found a constant ratio of gas phase 
H2CO/CH3OH between 0. 1 3 and 0.28 . 

In this paper, we aim to develop a model that explains these observed abundances and abundance ratios an d to 
study the impl i cation s for the conditions in regions where they are found. Following iTielens & Hagen d 19821) and 
Charnlev et al. (1992) and supported by nume rous laboratory studies on hydrogenation of CO (iFuchs et al] 120091: 



Hiraoka et al] l2002t IWatanabe & Kouchill2002l) . we assume that methanol and formaldehyde are formed from hy- 
drogenation of CO on grain surfaces. Observations have revealed that CO is present in three distinct components 
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in interstellar ices: non-polar (pure) CO, C O mixed with CO2 and CO as a trace in polar ice dChiar et alJll998c 
Pontoppidan et al.ll2008t iTielens et al.ll 1 99 lb . Here, we exclusively focus on methanol formation during the conden- 
sation/formation of the non-polar CO component. This component is thought to be formed when the cloud density 
increases and the gaseous atomic H/CO ratio decreases (ITielens & Hagenll 1982b . and thus occurs at a later stage or 
deeper in the cloud than the water ice formati on. This is consiste nt with the higher observed extinction threshold for 
CO ice formation compared with that of H2O dWhittet et al ] l200ll) . In the highest density regions (> few xlO 5 cm 3 ), 
the timescales for collisions of CO with the grains become so short that most of the gaseous CO is removed from the 
gas. This so-called 'catastrophic' CO freeze-out has been observed directly through infrared ice-mapping observations 



( Pontoppidanll2006T) and in directly through the lack of gas-phase CO (and other molecules) from the densest parts of 
the core, e.g., L1544, B68 ( Bergin et al. 2002 : Caselli et al.ll 19991 ). as well as the accompanying rise in H2D + and the 



deuterated molecules. In such high density regions, the bulk of the CO ice is in the pure form. For this pure CO-ice 
component, we can thus be assured that CO is the main reaction partner of accreted hydrogen and hence all the rele- 
vant reactions have been studied experimentally. For the other CO ice environments, competition with other reactions 
clearly plays an important role, in particular hydrogenation of atomic oxygen to water. We briefly comment on this 
competition in §5 but leave the formation of methanol and formaldehyde in these other interstellar ice environments to 
a future paper. 

The formation route of methanol ice by hydrogenation of solid CO studied in the laboratory consists of the follow- 
ing steps: 

CO -> HCO -> H 2 CO -> H 2 COH -> CH3OH, (1) 

where the reactions involving CO and H2CO possess small activation energy barriers. The experiments all show that 
for a range of temperatures and H-fluences (time x flux) that are comparable with fluences in the ISM, the H2CO rate 
of production is equal to or higher than the fo rmation rate of CH 3OH. This result appears to be in contradiction with 
the observational evidence. In a recent paper dFuchs et al. I l2009h . we showed by means of microscopic Monte Carlo 
simulations that high H-atom laboratory fluxes and low H-atom interstellar fluxes do not result in the same production 
rate of H2CO and CH3OH, even if the overall fluences are the same. The difference arises because most of the hydrogen 
atoms in t he high-flux regim e form molecular hydrogen whereas in the interstellar regime the hydrogenation reactions 



dominate (Fuc hs et al. 20091) . The crossover from F^CO-rich conditions to CH30H-rich conditions occurs therefore 



at much earlier t imes than expected from a direct translation of the laboratory environment. Moreover, the simulations 
presented in the Fu chs et al. I d2009h paper are for hydrogenation on pure CO ice, a situation which is unlikely in the 
ISM. Indeed, from the final interstellar yield of only 4 monolayers (ML) obtained after 2 x 10 5 yr, it was concluded that 
conversion of CO to CH3OH ice must occur simultaneously with the freeze-out and building up of the CO layer, since 
hydrogenation after freeze-out does not result in the high methanol abundances that are observed within reasonable 
times. Furthermore, although hydrogen resides in dense interstellar clouds predominantly in the form of H2, a fraction 
is present in atomic form together with gas phase CO, which suggests simultaneous hydrogenation and freeze-out. The 
present paper presents simulations of methanol formation under cold core conditions, where hydrogenation and freez- 
ing out of the CO both occur. These conditions also pertain to the outer envelopes of YSO's. Different temperatures 
and initial abundances are used and the results are compared with observational data on the H2CO/CH3OH gas and ice 
abundance ratios. The laboratory results will further be extrapolated to lower temperatures, which are more relevant to 
cold dense cores. 

The models presented in this paper differ from previous models on the formation of interstellar methanol ice in 
that they include detai ls of the surface structure of CO and utilise parameters obtained to reproduce experiments in 
the laboratory setting dFuchs et al.ll2.009r) . In this manner, a set of energy barriers for the different processes on the 
surface - diffusion, desorpti on and reaction - was obtained, which can now be used in a n interstellar environment. 
Unlik e the master equation dStantcheva et al.ll2002l) . rate equation ( Ruf fle & Herbst 2000), and macroscopic Monte 
Carlo dRuffle & Herbstll2000l) studies of methanol formation, the continuous-time, random-walk (CTRW) Monte Carlo 
technique accounts for the layering of the CO. This layering is crucial, because the constant addition of fresh CO tends 
to protect the underlying layers from hydrogenation and limits the time available fo r reaction. Moreov er, our current 
approach is different from previous CTRW-Monte Carlo simulations on this system ( Chang et alJl2007h . both because 
we use newly determined barriers, and because we account for the actual crystal stru cture of CO ice, instead of utilising 
a simple cubic structure. On the other hand, the previous study of lChang et alj d2007l) coupled the gas-phase chemistry 
with grain-surface chemistry and used a larger grain surface reaction network. 

The paper is organised in the following way. The next section describes the Monte Carlo method and introduces the 
input parameters of the model and their origin. Section[3]presents th e model results fo r interstellar conditions covering 
the same temperature regime as in the experiments (12.0-16.5 K) dFuchs et al]l2009l) . The focus of the discussion is 
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on how the abundances depend on various parameters, especially the gas-to-dust ratio, and on the ice structure and 
layering. Section [3~4l extends the model to lower temperatures of 8 and 10 K, which are more appropriate for cold 
dense cores. In Section |4]the model results are compared to a rate equation model using similar input parameters and 
conditions, an analytical steady-state model, other prior grain surface models, and observations. Section [5] discusses 
the effect of the limited chemical network on the results by introducing more reactions. 



2 Monte Carlo method 



The CTRW method has been described previously (s ee, e.g.,|C hang et al. 2005). A detailed overview of the Monte 
Carlo program and its input parameters can be found in lFuchs et al. (2009). In brief, the technique simulates a sequence 
of processes that can occur on a grain surface, which is modelled as a lattice with the number of lattice sites determined 
by the size of the grain and the site density for the adsorbate CO. The order of this sequence is determined by means of a 
random number generator in combination with the rates for the different processes. These processes include deposition 
onto the surface, hopping from one lattice site to a nearest neighbour, desorption of the surface species, and reactions 
between two species. We assume that the (first-order) rate coefficients (s _1 ) for hopping and desorption are thermally 
activated according to the formula 

/? x = vexp(-^) (2) 
where v is the pre-exponential factor, which is approximated by a constant value obtained fr om the transition state 



theory expression kT/h ~ 2 x 10 11 s 1 at low temperature (see Eq. (4.45) in KolasinskH 2002 ). and E x is the barrier 



(or simple endoergicity) for process X in K. These energies can be constrained by reproducing the laboratory results at 
different temperatures or obtained via ab-initio calculations. Unlike desorption and hopping, the rate coefficients for 
those chemical react ions that have an activation energy barrier also have a tunnelling component. In previous CTRW 
Monte Carlo studies ( Chang et al ,l2007l: E uppen & Herbstl2007h . we treated reactions with barriers solely by tunnelling 
through a square potential and ignored any small temperature dependence at low temperatures. The present paper u ses 
reaction rate coefficients that have a slight temperature dependence, which were obtained in Fuchs et al. ( 20091) by 
using the rate expression for thermal activation (Eq. [2]) and allowing the activation energy barrier for reaction to be 
temperature dependent. Table Q] lists the resulting barrier heights and reaction rate coefficients for the reactions H + 
CO and H + H2CO at temperatures from 12.0 K to 16.5 K. The barriers increase by ~ 20% from the lowest to the highest 
temperature. Since the intermediate products, HCO and H3CO, were not experimentally observed, it was concluded 
that the hydrogenation reactions of these species have a negligible barrier. In the program, these two reactions proceed 
without any barrier at all. After a reaction, the formed product gains some energy to rearrange and it is allowed to make 
a few diffusion steps. 



Table 1: Reaction rate coefficients, 7? leact , and activation energy barriers for CO + H and H2CO + H for different 
temperatures. 
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CO + H 




H 2 CO + H 






barrier 


^react 




barrier 


^react 
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(K) 


(s" 1 ) 




(K) 


(s" 1 ) 




12.0 


390 ± 30 


2x 10" 




415 + 30 


2x 10" 


-4 


13.5 


435 + 40 


2x 10" 


•3 


435 + 40 


2x 10" 


•3 


15.0 


480 + 50 


3 x 10" 


■3 


470 + 50 


5 x 10" 


■3 


16.5 


520 ± 60 


4x 10" 


■3 


490 ± 60 


2x 10" 


-2 



The desorption, or binding, energy for an adsorbate depends on an energy parameter E. To understand the role of 
this parameter, let us first consider a CO molecule either in the ice or adsorbin g onto the surface. The CO molecule 
is assumed to lie or stick in a configuration close to the a-CO structure dVegardlll930h . In this configuration, each CO 
molecule has up to 14 neighbours to which hopping can occur, four in the same layer, four one layer above and below 
and one two layers above and below. The binding energy £bmd is determined by the additive energy contributions of the 
occupied neighbouring sites. The contributions are 2E for the layers below and, if applicable, E for the neighbours in 
the same layer (lateral bonds) or in upper layers. An alternative treatment for sites below the particle is to add a single 
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Table 2: The energy parameter E for the different species. 



Species 


IT 

rl 


CO 








H 


3 


32 


H 2 


3 


33 


CO 


32 


63 


HCO 


160 


1600 


H 2 CO 


160 


1600 


H 3 CO 


160 


1600 


CH 3 OH 


160 


1600 



overall contribution for longer range interactions from the ice layer. In this way, the added long range contribution is 
roughly the same for all species of the same kind, regardless of their very local environment. However, if a very porous 
structure forms, the lower neighbouring sites are most likely not all occupied and the added long range contribution is 
less, reflecting the 'real' long range contribution. As an example, if a CO molecule lands on top of a CO layer in a 
multi-layered CO ice, its binding energy is l(2£'co-co)+4(2£ , co-co) = 10£co-co, while if it lies in a deeply embedded 
layer in the ice, its binding energy is the full l(2£ C o-co) + 4(2£ , C o-co) + 4(£co-co) + 4(£co-co) + l(£co-co) = 
19£co-co- From ab-initio calculations (Andersson, in prep) the parameter £co-co f° r CO is 63 K, so that the binding 
energy of a CO adsorbate onto CO ice is 630 K. Experimentally a desorption energy of CO from a CO surface was 
determined to be 855 + 25 K, which corresponds to 13-14£co-co- For atomic hydrogen, Eu-co was calculated to be 
32 K onto CO ice (A ndersson, in prep). For H 2 a value of £b,-co _ 33 K is used, slightly higher than for H as expected 



Bret) 

from calculations of iHornekaer et alj d2005l) for H 2 and lAl-Halabi & van Dishoeckl (12007b for H atoms. If lattice sites 



are occupied by heavier species than CO (mainly CH3OH) the total binding energy of an H or CO, and the structure 
are assumed not to change. If lattice sites are occupied by H atoms rather than CO, the lower value of £x-h is used for 
those relevant sites, so that the total binding energy for a specific CO is an expression containing both the parameters 
for CO and for H. For atomic H, the van der Waals interaction with another atomic H in an adjoining site is rather 
small, i.e., £h-h = £fi-co/10. The energy parameters Ex-co forH 2 CO and CH3OH and the intermediates, HCO and 
H3CO, are chosen such that these species neither desorb nor diffuse (see below) from the grain at the temperatures 
studied. All adopted values for E are summarised in Table [2] 

The same parameter E is used to help determine the hopping barrier from site i to an empty adjacent site, j. If sites 
i and j have the same binding energy, then the barrier for hopping is given by the equation 

E hop (i,j) = {E, (3) 

where £ is a parameter that determines how high the barrier is compared with the binding energy. Here £ = 8 is 
used, so that the b arrier height is 80% of the binding energy of a surface adsorbate (see above). As reported in 
Fuchs et alJ ([2009), the model lesults are relatively insensitive to the exact value of this parameter. In general, the 



effect of the diffusion rate was found to become larger with temperature, where faster diffusion results in less effective 
hydrogenation, since the residence time per site decreases and more H atoms will desorb. The value of ^ = 8 was 
chosen to result i n acceptable simulation times and to reproduce the desorption energy/hopping barrier ratio found by 
Katz et al. (1999) for H and H 2 on silicates and amorphous carbon. The same parameter ^ is used for all different 



species: since it is multiplied by the species-specific parameter E, the difference in hopping rates between species is 
automatically accounted for. For this choice of H and H 2 are mobile for T > 8 K and CO for T > 16.5 K. Diffusion 
from low to high binding sites is possible from even lower temperatures. 

If the binding sites i and j have different binding energies, we add a second term to the formula for the hopping 
barrier in order to have the forward and backward rates obey detailed balance. The overall expression then becomes 

Z7 /■ -v tip , AEbinddJ) ... 
£ho P 0, j)=%E+ , (4) 

where AEbfadtt j) is the difference in binding energy between the two sites. Diffusion into the structure of the ice is 
included as well. Minimum energy path calculations suggest that CO and H can swap positions, enabling an adsorbing 
H-atom to penetrate into the CO ice (Andersson, in prep). The barrier for this process strongly depends on the layer 
in which the H-atom is situated. In the simulations the barrier for this event is (350 + 5{z\ + z 2 )) K for an H-atom to 
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swap between layers zi and zi, where z, stands for the depth of the atom in ML. This swapping mechanism was only 
found to occur for the top few layers, where there is more space to accommodate the movement involved in the swap 
(Andersson, in prep). The minimum energy path calculations do not include thermal effects, which could promote 
this mechanism by lowering the barrier in the upper layers. In the current implementation, the layer dependence of 
the swapping barriers effectively limits the swapping process to the top few layers. Deeper into the ice, the swapping 
mechanism is not thermally accessible. Hydrogen diffusion throught the ice via interstitial sites was found not to 
reproduce laboratory results and is therefore not included. 

Thus, each species can in principle undergo up to 15 processes: thermal desorption and hopping to one of the up 
to 14 neighbouring sites, or, if sites are occupied, reaction with the species in that site. Reactions without activation 
energy {e.g., H + HCO — > H2CO) occur with 100% efficiency, whereas for reactions with barriers, we must simulate 
the competition between the efficiency of reaction and that of hopping out of the adjacent lattice sites. Some of these 
processes can be forbidden. Bulk species are, for instance, not allowed to desorb and not all species react with each 
other. The rate coefficients of all these processes are determined using Eq. [2] and the barriers as discussed above. In 
this way, diffusion and reaction are automatically treated competitively. 

In addition to reactions occurring via the Langmuir-Hinshelwood mechanism, which happen as a result of diffusion 
into adjacent sites, the program includes reactions occurring via the Eley-Rideal mechanism, in which a gas-phase 
species lands atop a surface species and reacts with it. The barriers for reactions are assumed to be independent of 
mechanism. In the diffusive case, however, if two reactants do not tunnel under or cross over the activation energy 
barrier, they have multiple additional opportunities until they diffuse away from one another. In the Eley-Rideal case, 
we allow only one opportunity. For example, an H atom landing directly atop a CO molecule gets one chance to 
overcome the reaction barrier and react. If the reaction does not proceed, the H atom can then hop further to find 
another reactant. 

Deposition onto a surface site occurs with a rate coefficient according to 



-ftdep - A . A , (5) 

4p 



where «a is the absolute gas abundance of species A, p is the surface site density, and va is the mean velocity of species 
A: 



VA = 



V 



(6) 



mnA 



with r gas the temperature of the gas and wa the mass of A. We use the standard value of p = 1 x 10 15 cirT 2 which agrees 
with the site density of the a-CO (110) crystal face. During the simulation H, CO and H2 are allowed to deposit on 
the surface. We use a molecular hydrogen abundance that is only slightly higher than the atomic hydrogen abundance 
to maintain a coverage of H2 on the surface within a reasonable simulation time. The presence of H2, which is also 
formed on the surface, has an effect on the penetration of the atoms and the diffusion. The CO gas-phase abundance 
is adjusted to account for freeze-out during the simulation. This freeze-out rate (loss rate of gaseous CO), which is the 
deposition rate onto a surface site multiplied by the number of sites per grain and the grain concentration, depends upon 
the dust-to-gas number ratio for fixed gas-phase density «e- This paper uses two values for this parameter: 1 x 10 _12 «h 
and 2 x 10~ 12 «h- The influence of the rate of depletion, i.e., the rate of disappearance of CO from the gas, is discussed 
in Section[32] Gas-phase chemistry is not included here. 



3 Results 



3.1 Ice abundances as functions of time 



In order to follow the build up of CO, H2CO, and CH3OH in the ice mantle, simulations at different temperatures and 
densities were carried out. All simulations commence with the sa me initial bare surfaces c onsisting of lattices with 
50 x 50 sites and some simple surface steps (similar to Surface c in lCuppen & Herb st (2005)) The choice of the initial 
surface (either silicates or amorphous carbon)) was found to be unimportant after the build-up of a few monolayers 
of ice. To start the simulation on a bare surface, the energy parameters used were chosen to be the same as for the 
CO surface. All results were converted to grains with a standard size of 0.1 yum. The simulation surface is a square 
flat surface. The spherical, continuous nature of a grains is mimicked by periodic boundary conditions. The si ze of 
the simulated surface has been found to be in the regime where no size effects are observed (IChang et al.l 2005): this 
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regime is associated with surface abundances of reactive species such as H larger than one per grain. For H2 formation, 
a size dependence is observed, since for smaller grains the average surface density of atomic hydrogen is much smaller 
than one (accretion limit), but because of the stronger sticking of CO and H2CO with respect to H atoms, the size 
dependence is probably only important for higher temperatures where CO and H2CO start to desorb. 

Figure Q] shows the evolution of the three dominant ice species as a function of time for n H = 1 x 10 4 cirT 3 , 
Wgrain = 2 x 10~ 12 n H and f our different tempe ratures: 12.0, 13.5, 15.0, and 16.5 K. These temperatures are the same as 
used in the experiments in Fuc hs et al. (2009). The abundances are plotted in terms of 10 15 molecules per cm 2 , which 
is equivalent to species per lattice site. The plotted curves are the combined results of several independent simulations 
using different initial seeds, e.g., initial settings of the random number generator. Simulations were added until the 
evolution did not significantly change. Initial gas phase abundance s of n(CO) = 1 x 10~ 4 «h and n(H) = 1 x 10~ 4 «h 
were assumed; these correspond to observed (CO) ( Lacy et all 1994 ) and calculated (H) fractional abundances for cold 
dense cores of a density of tin = 1 x 10 4 cirT 3 . The abundance of atomic H is mainly determined by the ration between 
the cosmic-ray destructio n rate of H2 and the form ation rate of H2 which leads to ~ 1 cm 3 independent of the density 
dDulev & Willia ms 1984). iGoldsmith & Lil d2005l) showed observationally that the H-atom abundance is higher than 
1 cirT 3 and we therefore use a higher H abundance of 10 cirT 3 for our high density case which will be discussed below. 
The gas-phase abundance of CO is allowed to deplete as CO accretes onto grains while the gas-phase abundance of H 
remains constant. The solid line indicates the combined ice build-up of the three species together; in 10 5 yr, at most 
10 monolayers of ice are produced at the value of nn used. From the graphs, a difference in this total build-up for 
the higher temperatures is immediately clear. Whereas the ice thickness grows linearly for 12.0 and 13.5 K, a clear 
non-linear behaviour as well as a drop in the total coverage can be observed at 15.0 and especially 16.5 K. This is due 
to the binding energy of CO. The residence time of an adsorbed CO molecule is 0.3 yr on a flat grain at 16.5 K. This 
time will increase once the CO molecules can stick together and form small islands, because the binding within the 
islands is stronger than for individual molecules. With the CO gas-phase concentration of 1 cm 3 , the arrival time of 
new CO molecules onto a grain is of similar order as the residence time on a bare surface, so that the CO molecules do 
not have the opportunity to meet and stick together to form small islands, which is an efficient mechanism for growth. 
Figure|2]shows similar time evolution curves for a density of «h = 1 x 10 5 cm 4 , leading to a ten times higher flux, and 
indeed here there is a greater ice build-up at 16.5 K. 

In general, the production of H2CO and CH3OH decreases with increasing temperature for both densities. Thi s 
dependence is in agreement with the laboratory experiments of the hydrogenation of a CO ice by Fuc hs et al. (2009). 
There the initial production rate was observed to be higher for low temperature whereas the final yield at the end of the 
experiments peaked at 15.0 K due to an increase of the penetration depth of H-atoms into the CO ice with temperature. 
In the co-deposition simulations presented in this paper, the penetration depth has less impact on the formation of H2CO 
and CH3OH since a fresh supply of CO is constantly deposited and the formation rate is therefore more comparable to 
the rate at the start of the experiments when a pure CO ice is exposed to atomic hydrogen. In this condition, it is the 
shortened grain lifetime of H atoms at the higher temperatures that reduces the rate of methanol formation from CO. 
The rate for reaction changes only minimally with increasing temperature. 

The pure CO build-up here peaks at 13.5 K for both densities as can be clearly seen inFigs.Q]and[2]by comparing the 
dash-dotted curves. This quantity is influenced by two effects with opposite temperature dependence: the desorption of 
CO molecules (see above), which increases with increasing temperature, and the conversion of CO into H2CO, which 
decreases with increasing temperature because of the lessened ability of H atoms to remain on grains long enough 
to react. Nevertheless, for most conditions, CH3OH dominates the ice layer at late times and exhibits a very steep 
formation curve. CO and H2CO on the other hand start with a high abundance at early times and increase much more 
slowly, even reaching a steady-state level on some occasions. The steady-state can be clearly seen for T = 12.0 and 
13.5 K for rin - 1 x 10 5 cm -3 and T = 16.5 K for riu = 1 x 10 4 cm 4 . For lower atomic H abundances, ~ 1 cm -3 , 
the H2CO and CH3OH formation is expected to occur at later times, since similar «(H)/n(CO) ratios will occur at later 
times. 



3.2 Dust-to-gas number ratio 

The dust-to-gas number ratio determines the rate of the gas phase depletion and the maximum ice thickness that can 
be achieved before depletion of gaseous CO, which is roughly 50 ML for n gra i n = 2 x 10 _12 «h and 100 ML for 
"grain = 1 X 10~ 12 «h as is discussed in the next section. 

Figure [3] shows the temporal evolution of CO, H2CO and CH3OH for 12.0 and 16.5 K with a reduced dust-to-gas 
ratio of 1 X 10~ 12 «h- Densities of «h = 10 4 and 10 s cm -3 are used. Panel (a) can be compared with the 12.0 K panel 
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Time [years] 



Figure 1: CO, H2CO, and CH3OH build-up as a function of time for a density of «h = 10 4 cm 3 , n gra i n = 2 x 10 12 «h 
and w(CO)/ra(H) = 1 at four different surface temperatures (a) 12.0, (b) 13.5, (c) 15.0 and (d) 16.5 K. The black line 
indicates the total ice thickness, i, e. the sum of the CO, H2CO, and CH3OH abundances. 



in Fig. [T^, panel (b) with the 16.5 K panel in Fig. Q]and panels (c) and (d) with the 12.0 and 16.5 K panels in Fig. [2] 
respectively. It is apparent that only Fig. [3J; is significantly different from its analogs. For the lower density cases 
(«h = 10 4 cm -3 ), substantial depletion of gaseous CO has not started yet after 2 x 10 5 yr, resulting in roughly the same 
gas-phase composition throughout the simulation for both dust-to-gas ratios. This is generally not true for the higher 
density models. At 16.5 K, however, the build-up of ice layers is hampered by the desorption of CO back into the gas 
phase, again resulting in very similar accretion rates for the two dust-to-gas cases. At lower temperatures, however, 
the accretion rate for CO is greater for the lower dust-to-gas case. The main difference between Figs. [2^ (12.0 K) and 
[3k is indeed that the levelling off to constant values for CO and H2CO occurs at later times and thicker ice layers for 

"grain = 1 X 10~ 12 H H - 

3.3 Ice structure 

In addition to the overall surface abundance as a function of time, the Monte Carlo approach allows us to obtain 
more detailed information concerning the layering of the ice. Fig.|4]plots the fraction per monolayer of the main ice 
components CO, H 2 CO and CH 3 OH at a time of 2 x 10 5 yr for 12.0 K (a, b, d)) and 15.0 K (c). The abscissa is 
the number of the monolayer starting from the bare surface as zero indicated by 'grain boundary' in the figure. The 
label 'gas phase boundary' indicates the top layer of the ice mantle, which faces the gas phase. Despite differences in 
dust-to-gas ratios and other parameters in the panels (see caption), the plot clearly shows that there is a gradient in the 
ice composition. At the onset of the CO freeze-out, a fraction of the layers remains in the form of CO and the H2CO 
is not fully hydrogenated to CH3OH. While the gas phase CO slowly depletes, the «(H)/n(CO) ratio becomes more 
favourable for the complete hydrogenation of CO. For example, in the 2 x 10 5 yr of the simulation, the gas phase CO 
abundance drops from 10 cirT 3 to 0.2 cm" 3 for 12.0 K, «h = 10 5 cm" 3 and n gla i n = 2 x 10~ 12 «h- The right panels 
show the ice composition for slower CO depletion (n gra i n = 1 x 10" 12 «h)- Panel (d) has additionally an altered initial 
abundance of n(H)/n(CO) = 0.5. Note that the horizontal scale for the right panels differs from that for the left panels. 
For n gla in = 1 x 10" 12 «h the change to pure CH3OH layers occurs over more layers than for n gra i n = 2 x 10" 12 «h- 
If panels (a) and (b) were plotted as a function of the fraction of the total ice thickness instead of the absolute ice 
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thickness, very similar graphs would be obtained. The final overall H2CO/CH3OH ratio after freeze-out is therefore 
independent of the exact dust-to-gas ratio. Comparison of panels (a) and (c) confirms that the rate and efficiency of the 
conversion of CO into methanol is determined by the temperature while comparison of panels (b) and (d) shows that 
the rate and efficiency are also determined by the n(H)/n(CO) ratio. A closer look at panel (d) compared with panel 
(b) shows that the lower n(H)/n(CO) ratio leads to larger H2CO and CO fractions for the lower monolayers and that 
the conversion to pure methanol layers occurs only at higher layers, which are formed at later times. 

Figure |5] paints a schematic picture of a grain mantle. At the onset of CO freeze-out, the flux of CO molecules 
accreting onto the grain is large, larger than the part of H atom flux that is available for hydrogenation of CO, and, 
as a consequence, the mantle will be rich in CO with H2CO and CH3OH as minor components. As the n(H)/«(CO) 
ratio in the gas phase increases, more and more CO and H2CO is converted into CH3OH and the outer layers of the 
mantle become more CH3OH rich. Note that the astronomically observed solid abundances are sums over the entire 
grain mantle, although different ice components can be distinguished through the line profiles. 

Actual simulated grain mantle cross sections for individual runs are plotted in Fig. [6] Here unoccupied sites are 
indicated in black, light gray on the bottom of the figures represents the bare grain, while H and H2 are represented by 
white. All other mantle species, CO, HCO, H2CO, H3CO, and CH3OH, have different Gray scale levels according to 
their degree of hydrogenation; i.e., CH3OH is darkest, whereas CO is represented by the lightest Gray. The four panels 
correspond to the resulting icy mantle after 2 x 10 5 yr under the same conditions as in Fig. [2] The grain temperatures 
are 12.0, 13.5, 15.0 and 16.5 K from top to bottom. The zigzagging pattern at the grain mantle surface reflects the 
zigzagging structure of a-CO. The panels show clearly the gradient in CO, H2CO and CH3OH across the grain with 
the top layers containing the more saturated species. Furthermore, it shows that the lower temperature grains contain 
more CH3OH. The grain mantle at 16.5 K is much thinner due to the desorption of CO at this temperature. Finally, some 
small pores can be observed in the mantles; these appear to form during hydrogenation where two species recombine 
to one, which takes up less space. It appears that at late times, when the flux is low, either CO or CH3OH has some 
time to rearrange and fill most of these vacancies. 
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Figure 3: CO, H 2 CO, and CH 3 OH build-up as a function of time for n grain = 1 x 10 12 n H and «(CO)/n(H) = 1. (a) 
T = 12.0 K and ra H = 10 4 cm" 3 , (b) T = 16.5 K and n H = 10 4 cm 4 ,°(c) T = 12.0 K and n H = 10 5 cm 4 , and (d) 
T = 16.5 K and n H = 10 s cm -3 . 
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Figure 4: Fractional abundance of the three main ice components as a function of monolayer after 2 x 10 s yr at «h = 
10 5 cm -3 for (a) 12.0 K and « glain = 2xlO" I2 « H , (b) 12.0 K and n grai „ = lxlO" 12 n H , (c) 15.0 K and « grain = 2xlO" 12 « H , 
and (d) 12.0 K and n gla i n = 1 X 10 _12 «h with a lower initial n(H)/n(CO) gas phase abundance ratio of 0.5. 
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Figure 5: Schematic picture of the growth of the ice mantle during CO freeze-out. For coloured figures, see the online 
version. 



3.4 Methanol formation at 8 and 10 K 

The results presented above all concern the formation of methanol in the temperature regime between 12.0 and 16.5 K. 
These temperatures are more representative for the outer regions of protostellar envelopes. In cold dense cores, the 
temperature is most likely lower. Figure [7] presents simulation results for 8 and 10 K. The model is extended to 
temperatures outside the regime studied in the laboratory, by assuming that the reaction rates (/? rea ct in Table [TJ do not 
change with temperature, since they are dominated by tunnelling, whereas all other processes are thermally activated. 
The results in Fig. [7] clearly show that formaldehyde and methanol are still efficiently formed at these temperatures. 
Simulations with much lower reaction rates for H+CO and H+H2CO, thermally activated using barriers of 390 and 
415 K, respectively, still result in the formation of both H2CO and CH3OH, although with much lower abundances 
as indicated by the thin lines in Fig. [7] The abundance of formaldehyde is much larger in this case, since it is not 
efficiently converted into methanol. However, the rates at higher temperature suggest tunneling to be important and the 
abundances are therefore probably better represented by the thick curves. 



4 Comparison with other models and observations 

This section compares the results of our Monte Carlo simulations with different models and with observations of H2CO 
and CH3OH abundances. In particular, our result s are compared with a similar model with similar input parameters 
using a standard rate equation technique (see, e.g., Ruffle & Herbsdl2000l) and with different models of various degrees 
of chemical and astrophysical compl exity that are reported in the literature. Finally, the simulations are compared 
with a steady-state model proposed bv lCharnlev et al.l ( 1 1997b that we have adjusted to take the changing CO gas-phase 
abundance into account. 



4.1 Comparison with rate equations 

Figure [8] plots the results for a rate equation model (see, e.g., Ruffle & Herbst 2000l) with the same chemical and 
physical processes as used in our Monte Carlo approach. Besides the difference in mathematical techniques, the rate 
equation method uses single hopping and desorption rates rather than rates that depend on the local structure. We 
expect the diffusion to be dominated by hopping between sites of the same type which results in a diffusion energy 
of £hop = 8£ (see Eq. |4|i and the desorption dominated by strong binding sites (three horizontal neighbours), since 
diffusion will allow the particles to move to these sites where they remain attached. The competition between reactions 
with activation energy barriers and hopping and desorption is treated by dividing the rate coe fficient for reaction by the 
total rate including diffusion and desorption, as discussed in detail by lHerbst & Mill ar (2008). The panels in Fig.[8]can 
be compared with those in Fig. [2] which are determined with the Monte Carlo approach. 

In comparison with the Monte Carlo results, two trends become immediately apparent in the rate equation results: 
(1) there is a clear drop in the CO and H2CO abundances at late times and (2) the difference in the initial CH3OH 
formation (< 10 3 yr) increases with surface temperature. The drop in the surface abundance of CO and H2CO at late 
times is due to the mean field character of the rate equation technique. It treats all molecules of the same species in 
the same way, regardless of their position in the ice layers. CO ice that resides at the lower layers of the ice can in the 
standard implementation of the rate equation technique still react with impinging hydrogen atoms, whereas in reality 
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this reaction will not proceed unless H atoms can penetrate into the porous ice. Figure|4]shows that indeed most of the 
CO and H2CO is buried deep into the ice. In the Monte Carlo simulations, the hydrogenation reactions will therefore 
be hampered by the deficiency of gas phase CO at late times whereas the rate equations continue hydrogenating the CO 
and H2CO that have been formed at early times. In general, the effect of a changing gas phase composition is limited 
to the top layers of the ice. For this reason, laye ring should be taken into a ccount when modelling the grain surface 
chemistry. The three-phase model introduced bv lHasegawa & Herbstl d 1993b is a first attempt to implement this effect 
into a rate equation model. 

The second effect is probably due to the difference in the treatment of diffusion and desorption between the two 
methods. The agreement in the total ice thickness indeed becomes less if lower desorption barriers are used. Due to 
the implementation of the competition for reaction in the rate equation approach, the choice of the diffusion barrier 
has very little effect on the final results. This can be understood by realising that faster diffusion will lead to more 
encounters between diffusers and the more stationary reacting species but that the residence time in the vicinity of 
the stationary reactant is reduced. The site dependent rates in the Monte Carlo code result in longer residence times 
at some sites and shorter residence times at other sites, favo uring reactions, especially at higher temperatures, as first 



discovered for the formation of H2 jCuppen & Herbstll2005 ). In addition to the two effects, it can also been seen that 



with the rate equation approach, the ice cannot even develop to 1 monolayer at 16.5 K, presumably because the CO 
growth mechanism of forming islands is not accounted for. 



4.2 Comparison with other surface models 

A variety of different models concerning the formation of formaldehyde and methanol via surface reactions from 
CO have been reported in the literature. Figure [9] (top) shows fractional methanol and formaldehyde abundances with 
respect to «h as functions of the product of time and density («h) for a number of these models. This product is roughly 
proportional to the total fluence, the number of reaction species that reaches the grain surface, and can therefore serve 
as a measure to help compare models with different densities. The solid lines indicate the methanol and the dashed lines 
the formaldehyde abundances obtained in the present paper (Figs. [1] [2] and[7j. Only the data at 8.0, 12.0 and 16.5 K 
are plotted since the 12.0 and 16.5 K data represent the extreme values and 8.0 K is closed to cold core conditions. The 
«h = 10 4 and 10 5 crrT 3 data overlap nicely indicating that the formation rate of methanol and that of formaldehyde are 
directly dependent on the fluence and that the flux difference of one order of magnitude does not introduce additional 
scaling effects. 

The symbols represent the fractional abundances of both species obtained from a selection of different studies. 
Open symbols represent H2CO whereas filled symbols represent CH3OH. Given the large differences in parameters, 
and assorted met hods of calculation amo ng the various studies, large variance in results can be expected. The diamonds 
show results by Ruf fle & Herbst (2000), who used a large reaction network of both gas phase reactions and grain 
surface reactions. Both chemistries were treated using rate equations for four different scenarios defined by slow 
and fast diff usion rates combined w ith atomic and molecular initial conditions. The circles represent master equation 
results from Stantc heva et al. (2002). Here for a wide range of densities (10 3 , 10 4 , 10 5 cm" 3 ) the gas-grain chemistry 
for a limited set of surface reactions is obtained. Five additional reactions with respect to the surface network of the 
present paper were included, leading to the forma t ion of H2O, CO2 and O2. No gas phase chemistry was considered. 
The two triangles are obtained from Garro dTt alJ (2006), who followed both the gas and grain chemistry during the 
collapse and heat-up phase of hot cores using a rate equation approach. The points used here are in the early times 
of the collapse when t he density is still close to the initial density. Again a ful l network was used, comparable to 
Ruffle & Herbst (2000), but with intermediate diffusion rates. In all three models dGarrod et al 2006: Ruffle & Herbstl 
2000; IStantcheva et alJl2002l) the barrier crossing of surface reactions with activation energy was treated simply by 
multiplying the meeting rates of reactive species by the probability of crossing the activation energy barrier. Unlike 
our rate equation approa ch (see above), surface reactions with barriers were not treated to be in competition with other 
processes like diffusion dHerbst & Millarj|2008l) . This neglect magnifies the differences obtained between low and high 
diffusion rates. Furthermore, these three models do not capture the layered structure of the ice mantle but allow all 
species to react with each other, regardless of their position. This assumption can become a problem once multiple 
layers start to build up, as shown in the previo us section. 



The squares in Fig. |9] represent data fromlChang et alJ (120071) obtained by a similar Monte Carlo method to that 
used in the present paper dChang et alj|2005l) . which includes both the layering effect as well as the competition of 
the reaction barriers. The surf ace chemistry is coupled to the gas phase chemistry and the surface chemistry network 
is very similar to one used by IStantcheva et al.l ([2002). The square symbols represent results for «h = 2 x 10 4 crrT 3 
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and a surface temperature of 10 and 15 K. The main difference between the method by IChang et al.1 d2007l) and the 
present simulations is that our Monte Carlo algorithm is optimised to reproduce laboratory experiments and it includes 
important features of the CO + H system that were revealed through this optimisation. 

In general, Fig.|9]shows a large spread which reflects the different levels of complexity among the different methods. 
Most of the obtained abundances of H2CO and CH3OH lie between our 12.0 and 16.5 K results and follow roughly 
the same trend that H2CO dominates at early times and CH3OH at late times. It appears that the slow diffusion 
results underestimate the formation of both molecules and the CH3OH/H2CO ratio, as compared to our model that is 
parameterized by experiments. Our 8.0 K data lie on the high end of the results, probably for two reasons. Our n(CO) 
is relatively high throug hout the entire simulation, although it should be comparable to the intermediate abundance of 
Stantcheva et al.l d2002l) . and our limited network might overproduce H2CO and CH3OH instead of forming H2O and 
CO2 for instance. Our simulations are however geared towards conditions where most of the atomic oxygen is locked 
up in CO, and H and CO are indeed the most important reactants on the surface. 

The bottom panel of Fig. |9]shows the H2CO/CH3OH ratio for the low temperature models (< 20 K) presented in the 
top panel. It is immediately evident from this plot that there is a spread in this quantity for the different models, even 
within the same paper. The difference in the H2CO/CH3OH ratio for our 8.0, 12.0 and 16.5 K is much smaller compared 
with the spread found in other models. This spread is due to different assumptions in the modelling method and 
differences in methodology. In a few instances, the spread is due to different physical conditions such as temperature 
(as in our data) and density. In general, the largest differences with our results come from models with a low diffusion 
rate that do not consider the competition between reaction and diffusion. 



4.3 Comparison with analytical model 

Charnl ev et al ] (1 19971) derived an analytical steady-state rate equation approach to predict the CO/CH3OH and H2CO/CH3OH 



surface abundance ratios as functions of an, the relative H-to-CO flux ratio, and <pn — Pco/Ph,co , the ratio of the 
probabilities for single CO and H2CO molecules to react with an H atom; i.e., the ratio of the rate coefficients. The 
exact derivation of these express ions is given in Appen dix A. Each value for an and 4>n results in a unique combination 
of both abundance ratios. When lCharnlev et al.l (1 19971) was published, (pn was unknown. Now, we can use TableQ]to 
obtain this quantity for different temperatures, so as to compare this steady state rate equation approach with the more 
realistic Monte Carlo simulations. 

The top panel of Fig. [TO] plots the n(CO)/«(CH30H) abundance ratio versus the n(H2CO)/n(CH30H) abundance 
for the Monte Carlo simulations (open symbols except for 16.5 K) and for the steady-state model (lines). For each 
temperature, (pn is indicated as obtained from Table Q] The Monte Carlo curves follow the abundance ratios in time; 
i.e., they start from a relatively high CO gas-phase abundance which is partially converted into formaldehyde and 
methanol, so that the general direction is right to left. Likewise, the general direction is from top to bottom since 
formaldehyde is generally converted into methanol. The steady-state model by ICharnley et al.l (119971) plots curves for 
varying an- Since an increases gradually during the course of our simulations, both approaches can be compared in 
a relative straightforward way although the simulations have also other sources of time dependence. The time it takes 
for the top layer of the grain mantle to reach steady state is short as compared to the change in an- Figure ITOT top) 
includes the steady-state mode l lines for values of 4>n that correspond to the values in the simulations. We would like 
to emphasise that the lines by ICharnley et alJ (1 19971) represent 0co/#ch 3 oh and #h^co/#ch,oh, where 9x indicates the 
steady-state coverage of species X in the top layer, whereas the abundances in the Monte Carlo simulations are for the 
entire grain mantle. At first sight, very good agreement is obtained between the simulation results and the analytical 
model, especially considering the requirement for steady state in the latter, which is not fulfilled in the Monte Carlo 
simulations. However, it must be remembered that in the absence of a third dimension to the plot, an is a hidden 
parameter, and the agreement for an between the two methods for the same points on the figure is not good. Indeed, 
the same results on Figure [10] are obtained for a 5-20 times lower H/CO flux in the steady-state model as compared 
with the simulations. In the simulations, a signifi cant portion of the H atoms are 'lost' to H2 formation and desorption. 
Both processes are not included in the model bv lCharnlev et alJ (11997b . 

The smaller symbols in the right top corner are the result of simulations with a different starting H abundance, 
chosen such that the flux of CO molecules is 10 times higher than the H atom flux, e.g., a//is 0.1 at the start of the 
simulation. The absolute abundances of H and CO have little influenc e on the obtained ratios. These independent 
simulations continue nicely on the analytical curve. As iKeand d200lh shows, the introduction of more competing 
hydrogenation reactions leads to different curves. Section|5]addresses this point in more detail. 

One purpose behind a plot like Fig.[l0]is to visualise whether H atoms are more likely to react with CO or H2CO. 
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The laboratory results indicate that for low temperature (12.0 K) CO + H is dominant whereas H2CO hydrogenation 
becomes the preferred channel for higher temperatures (16.5 K). The plot shows that, even though the ore agreement is 
not very good, the agreement for <pn is very nice. The curves in Fig. lTOh op) are for one value of 0h and varying an, 
which can be compared with the simulations at a single temperature. The steady-state model and Monte Carlo simula- 
tion results for the same value of <p H trace the same unique lines in the n(CO)/n(CH30H) versus n(H2CO)/n(CH30H) 
plot. Since H decreases with temperature, plotting observational H2CO/CH3OH versus CO/CH3OH ratios could give 
a good temperature and evolutionary indication between sources. 

The main conclusions that can be drawn from these comparisons are that the surface chemistry of the top layer can 
at each point in time be very well approximated by quasi-steady state conditions, since the time to reach a quasi-steady 
state is smaller than the change in flux, and that the majority of the hydrogen atoms do not react but leave the surface 
by desorptio n. This means that the conversion rate of CO and H2CO into more saturated species is overestimated by 
the model of ICharnlev et al.l (11997b . Table[3]lists the steady-state conversion fractions of CO into H2CO and CH3OH 
as a function of temperature and H/CO gas phase abundance ratio obtained from the Monte Carlo simulations. The 
calculated fractions include the effect of desorption of atomic hydrogen and the competition of the CO hydrogenation 
reactions with the formation of H2, which is not included in the analytical model. 

These conversion fractions clearly show that CO is very efficiently converted into H2CO and CH3OH. Naturally, 
the exact numbers depend on many assumptions regarding for instance the values of E, £ and the reaction barriers. The 
parameter E is central in determining the desorption temperature of the species and an increase in E will therefore result 
in an appreciable H-atom abundance at higher temperature increasing the hydrogenation regime with a few degrees. 
The diffusion parameter on the other hand, will have a stronger effect at low temperature. Finally, a change in the 
reaction barriers will lead to a slightly different H2CO/CH3OH ratio. 



4.4 Comparison with observations 

A comparison between the simulation results presented in this paper and observations is most straightforward with IR 
ice data. As long as the ice is not processed too severely, either by UV photons, which can break down formalde- 
hyde and methanol, or by heating, which can sublimate CO, the IR data represents the most direct comparison with 
our simulations (see Section [5] below). Ideal conditions are most likely found in cold dense cores or the outer cold 
envelopes of YSO's. Since it is difficult to use absolute column densities for the comparison, the same ratios are used 
as in Fig. [lOl CO is usu ally present on grains in different mixtures: pure CO, polar CO and CO mixed with CO2 
dPontoppidan et al.l l2008). Laboratory spectra of CO mixed with either H2O or CH3OH show similar changes in th e 
4.7 /xm CO feature, broadening and redshifting ( Bisschopll2007 ; Bottinelli et alJlinpreparation ; Bouwman et al.ll2007h . 
The polar CO component at 2139 cm -1 , normally ascribed to CO in a water-rich ice, can therefore also be due to a 
methanol-rich mixture. A CO:CH30H= 1 : 1 mixture is already sufficient to shift the band to the observed range. Thus, 
to obtain the observed CO/CH3OH ratio for comparison with our models, we use both the pure CO component and 
the polar CO component to account for the maximum amount of CO mixed with the formed CH3OH ice. The trian- 
gles in the bottom panel of Fig. [TOl represent ratios obtained f rom ISO data for the YSO's AFGL 989, AFGL 2136, 
W33A, AFGL 7009S, and NGC 7538 IRS9 dGibb et aljbool and VLT-ISAAC data for the Class source Serpens 
SMM 4 (Pontoppid an et alJl2004h. A rrows indicate upper limits. As with the earlier data, Spitzer spectra of sources 
near SMM4 by lBoogert et al.l (120081) resulted only in upper limits for H2CO. They f ound forma l dehy de to contribute 
perhaps 10-35% of their CI component, which cannot be assigned to a single species. iGibb et al.l ( 120041) us ed a feature 



equiva lent to the CI component to obtain their report ed H2CO abundanc es. Both the H2CO abundances by Gibbetal 



(2004) and the H2CO part of the CI components bv lBoogert et al.l (|2008) vary minimally with respect to water ice and 
are typically 6 %. The range in H2CO/CH3OH abundances is therefore mainly due to methanol. The abundance ratios 
for the six sources nicely overlap with the simulation data (solid lines). The good agreement with the CO/CH3OH ratio 
could however be deceiving if part of the CO ice has already desorbed. 

Solid methanol abundances ar e found to vary substantially with respect to water ice, ranging fr om upper limits of 
a few percent to more than 30% dBoogert et al.ll2008t iDartois et al.lll999t iPontoppidan et al.ll2003l) . In Fig. \W\ these 
percentages are indicated for the six sources compared here, spanning a range of more than an order of magnitude, 
assuming that, as discussed above, H2CO ice does not vary much in abundance. A reason for this large range could 
lie in the differing evolutionary stages of the sources. If most of the water ice is formed first before the catastrophic 
freeze-out of CO, from which methanol is formed, the methanol over water ratio will increase in time. The temporal 
evolution of the simulated abundance ratios proceeds from the top-right corner to the bottom-left corner of the figure, 
and the increase in methanol abundances for the six sources appears to follow this trend, suggesting that indeed the 
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Table 3: Conversion fraction of CO into H2CO and CH3OH as a function of T and H/CO 
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large spread in methanol observations can simply be due to a difference in evolutionary stage of the outer envelope. 
Unfortunately, solid H2CO has only been detected in a handful of sources and the same holds for clear upper limits for 
this molecule. 

An alternative observational test could be to perform CO, H2CO, and CH3OH gas phase observations in very cold 
and dense cores, where all three spec ies are frozen-out onto the grains. Since most molecules have similar non-thermal 
desorption rates dOberg et al. 2009allbl). o ne would expect the trace amounts in the gas phase to be representative of the 
ice layer composition ( Oberg et alj|2009h . The problem with this technique is that not all observed gas-phase CO may 
result from CO ice evaporation and that gas-phase reactions can contribute to H2CO as well. As shown in Figure 12 of 
Fuchs et al. j2009h . the current model can reproduce the observed H2CO/CH3OH ratios in high-mass hot cores where 



the majority of the observed gas-phase molecules are likely evaporated ice species. 

By ignoring the formation of water ice, our model predicts that CH3OH is mainly present on the grain mantles 
in the pure form or mixed with CO, and that it is not in a water-rich phase. Unfortunately, the ice composition has 
very little effect on the 9.75 jum band profile of CH3OH, which is usually used to determine its abund ance, neither in 
peak shape nor position. The feat ure only becomes red-shifted when water ice is dominant, >90% ( Bot tinelli et al.1 



in preparationy Skinner et a i1ll992l). Similarly, the 3.54 ^m CH3OH feat ure can be used to constrain the CH3OH ice 
"till 



environment ((Dartois et alJll999l:|Pontoppidan et alj2003l:lThi et alj|2006l) . All of these studies generally conclude that 
at least a fraction of the CH3OH ice is in a water-poor, CITjOH-rich environment but additional laboratory data on 
CH3OH mixtures with CO are needed to quantify this. 



5 Impact of competing reactions 

The Monte Carlo simulations reported in the previous sections all involve hydrogenation of pure CO ic e. The reason 
for this is two fold. First, one a im of this study is to model CO freeze-out in the centre of cold cloud cores dPontoppidanl 



2006; Pontoppi dan et al.ll2008l) where most of the gas is in molecular form and H2 and CO are the dominant species. 



In these centres, most of the elemental oxygen is in the form of CO or frozen out into grains in form of H2O. 

To study the effect of competing reactions, however, simulations that include oxy gen atoms have been performed 



as well. These simulations use a reaction network that is similar to that reported by IChang et alJ (120071) . leading to 



the formation of O2, H2O, and CO2. In addition, H2CO and H3CO can be destroyed by reaction with OH, leading 
to the formation of water (see Table HJ. All new species are assumed to bind very strongly and have only minimal 
diffusion. Eley-Rideal reactions are again allowed. The gas-phase abundances of CO and O are both assumed to be 
5 x 10- 5 n(H). Figure [TO] contains results for the simulation at 12.0 K using filled circles. From this figure, two things 
are immediately clear: the time dependence of the CO/CH3OH and H2CO/CH3OH ratios is minimal in the simu lation 
and th e filled circles overlap nicely with the analogous simulations that do not include the extra reactions. iKeand 



(120011) ^resents a similar graph obtained with a full gas-grain network using macroscopic Monte Carlo simulations. 
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Table 4: Surface reactions in the H, O and CO system 



Reactants Products E a (K) 

H + H ^ W 2 

H + O -> OH 

H + OH H z O 

H + CO HCO TableQ] 

H + HCO -> H 2 CO 

H + H2CO H3CO TableQ] 

H + H3CO -> CH3OH 

O+O -> 2 

O + OH -> 2 + H 

CO + OH -> CO2 + H 176 

OH + H 2 CO -> HCO + H 2 

OH + H 3 CO -> H 2 CO + H 2 



Their curves obtained for the full network are less steep than for the analytical expressions that are based on a limited 
network, since a considerable fraction of the hydrogen atoms reacts with other grain species. This is in contrast with 
our simulations that include H2O and CO2 as competing mechanisms. The origin of this discrepancy is unclear. 

Another type of competing reaction would be the destruction of methanol by photodissociation. The photofrag- 
ments could then react to synthesise more complex molecules. In dense cloud conditions , photodissociation mainly 



occur s via cosmic ray induced photons, with a typical flux of 5 x lCr photons cm s 



1992). 



Using a photodissociation cross section of 1.6 x 10 18 cm 2 dGerakines et al 



1996; 



( Cecchi-Pestelli ni & Aiellol 
Ober g et alJl2009h . 5 % of 



the methanol molecules has been photodissociated in 2 x 10 5 years. Most of these molecules are dissociated to simpler 
species like CO, HCO, CH3, CH3O, and H2CO and can be hydrogenated to methanol again. Gas-phase H atoms can 
hydrogenate the radicals in top layers of the ice, whereas H atoms that are produced through photodissociation can 
react with the fragments deep in the ice. We therefore expect only a minor change in the grain abundances due to 
photodissociation as long as the grains remain cold. Higher tempera tures allow the radicals to diffuse rap idly enough 
to form large molecules such as methyl formate and dimethyl ether (iGarrod et al feoOSHOberg et al.ll2009h . 



6 Conclusions 

The surface formation of CH3OH and H2CO from precursor CO has been simulated using the continuous-time, random- 
walk Monte Carlo method. The formation of both species was found to be very efficient under certain conditions and 
to depend mainly on grain temperature and the gas phase abundance ratio of H and CO. During the freeze-out of CO 
onto the grain, this ratio changes, favouring the more complete hydrogenation of CO to CH3OH. The more unsaturated 
species remain locked in the lower layers of the ice mantle. Due to the layering of the ice, changes in the gas phase 
abundance can only affect the top layers of the grain, which is important to take into account when modelling grain 
surface chemistry. The detailed Monte Carlo method used can be compared wit h a variety of other approaches. Of 
these, perhaps the most successful is a quasi-steady-state rate equation approach ( Charnlev et alj[l997 ) that focuses on 
the outermost layer of the grain only. 

The model results can be compared with observations through plots of the CO/CH3OH and H2CO/CH3OH abun- 
dance ratios. A very good agreement is obtained for the outer envelopes of a number of YSO's, with temperatures in 
the range 12.0-16.5 K. Moreover, the comparison allows us to trace the temporal evolution between different sources. 
Our results suggest that the difference in CH3OH abundances with respect to water are mostly due to differences in 
temporal evolution where the younger sources have not yet had sufficient time to build up much methanol. 
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A CO hydrogenation in steady state conditions 



This appendix presents a derivation of the steady state model proposed by ICharnlev et al. ( 1997 ). When the flux 



of hydrogen atoms to the surface is lower than the flux of CO molecules, the assumption can be made that every 
hydrogenation atom either reacts with CO or with H2CO. The probabilities of reaction are denoted by xco and;^H,co, 
respectively. They are determined by the likelihood of barrier crossing once an H atom is in the vicinity of the CO or 
H2CO, Pco and Pn 2 co, and the fractional coverage of these species over the outermost surface layer, 6qo and #h 2 co 
according to the equations 

co fcoflco ^. 

^CO#CO + ^H 2 C0#H 2 C0 

and 

Ph 2 co Oh, co 

XH 2 CO - — • (8) 

^CO^CO + ^H 2 CO#H 2 CO 



The change in CO coverage with time is then 

dflco 
dt 



= /co(1-0co)-2/h^co. (9) 



Here the first term accounts for the increase by the incoming CO flux, /co where the -6co accounts for covering of 
surface CO by new CO and the second, last term is due to reaction of CO with H. The factor two accounts for the 
two successive hydrogenation reactions to form H2CO. In a similar fashion, rate equations for the surface coverage of 
H2CO and CH3OH can be determined to be 

d#H,CO 

— — — = 2/hOtco - ^h 2 co) - /co0h 2 co (10) 

and 



dt 

From Eqs.[10]and[TT] #h 2 co/0ch 3 oh is given by 



d#CH 3 OH 

- 2/hVh 2 co - /co0ch 3 oh- (11) 



#H 2 CO , #CO 



#CH 3 OH #H 2 CO 

with 

, _ Pco 



Pi 



H 2 CO 



and furthermore 

#co 6h 2 co 6co 



#CH 3 OH #CH,OH #H 2 CO 

From Eqs.|9]and[T0]and using Eqs. [7][8j we obtain that 

Geo 1 + 20h/«h - (pn + V(l + 20H/O-H - <Ph) 2 + Hu/an 



&n 2 co 2</>h 

where 



(12) 

(13) 
(14) 

(15) 



an = (16) 
/co 
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Figure 6: Cross section of the ice mantle after 2 x 10 5 yr at a density of «h = 10 5 cm -3 for 12.0, 13.5, 15.0 and 
16.5 K (top to bottom). Unoccupied sites are indicated in black, grain by light grey, intermediates by cyan, and H and 
H2 by white. CO, HCO, H2CO, H3CO, and CH3OH are coloured by increasing darkness according to their degree of 
hydrogenation (see Fig |5j. For coloured figures, see the online version. 
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Figure 7: CO, H2CO, and CH3OH build-up as a function of time for a density of «h = 10 5 cm 4 , « gra i n = 1 X 10 -12 «h 
and «(CO)/«(H) = 1 at two different surface temperatures (a) 8.0, and (b) 10.0 K. The black line indicates the total 
ice thickness, i. e. the sum of the CO, H2CO, and CH3OH abundances. The reaction rates are assumed to be either 
governed by tunneling (thick lines) or purely thermal (thin lines). 




Figure 8: Similar to Figure[2]but using rate equations instead of Monte Carlo simulations. 
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Figure 9: Comparison of the present simulation data with values obtained by different models (see legend) as a function 
of time x density, which is a measure of fiuence. The top panel displays the H2CO and CH3OH surface abundances 
relative to «h, the bottom panel H2CO/CH3OH ratios. Lines represent the present model data for 8.0 (black), 12.0 (red) 
and 16.5 K (green) and «h = 10 4 cirT 3 and 10 5 cm 4 (overlapping). H2CO is in dashed lines and CH3OH in solid lines. 
In the top panel, the open symbols represent H2CO, the filled symbols CH3OH. 
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Figure 10: CO/CH3OH versus H2CO/CH3OH ice abundance. The top panel compares results as obtained by the Monte 
Carlo simulations (open symbols) and Eqs. [T2lfT4l (solid lines). The filled circles represent a simulation at 12.0 K that 
includes an O-atom flux. The bottom panel indicates the Monte Carlo results by solid lines. The triangles indicate 
ratios from ice observations where the numbers refer to the methanol content in the ice with respect to H2O (see text). 
The quantity 0h represents the ratio between the H+CO and H+H2CO reaction rates. 
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